#import gff for anno-TSS
features <- rtracklayer::import("Input/20210217_syne_onlyUnique_withFeat.gff3")
TUs <- rtracklayer::import("Input/Kopf_4091_TUs_combined.gff3")

1 Create .grps

Coverage data was created using Filtered .bam files with multireads, which were processed by bedtools genome coverage. Files were re-named in the respective history and downloaded as one .zip. In a next step, separate files were combined into one large table using Python code (joinFiles_PSS.py, joinFiles_TSS.py, joinFiles_transcript.py).

Normalization factors are size factors determined DESeq2 using separate .Rmd files (PSS-TSS Analysis.Rmd, TranscriptAnalysis.Rmd).

transcript_factors <- read.table("Input/transcript_sizeFactors.csv", sep=",", header=TRUE, row.names=1)
TSS_factors <- read.table("Input/TSS_sizeFactors.csv", sep=",", header=TRUE, row.names=1)
PSS_factors <- read.table("Input/PSS_sizeFactors.csv", sep=",", header=TRUE, row.names=1)
row.names(TSS_factors) <- c("dWT1", "dWT2", "dWT3", "WT1", "WT2", "WT3", "TV1", "TV2")
row.names(PSS_factors) <- c("dWT1", "dWT2", "dWT3", "WT1", "WT2", "WT3", "TV1", "TV2")
coldata_general <- read.table("Input/colData.csv", sep=",", header=TRUE, row.names=1)

1.1 PSS data

PSS_coverage <- read.table("Input/multireads_PSS_5ends_combined_5sensing.txt", sep="\t", header=TRUE, row.names=1)
PSS_fact <- PSS_factors$factor

if(FALSE){
ddsMat_PSS <- DESeqDataSetFromMatrix(countData = PSS_coverage,
                                 colData = coldata_general,
                                 design = ~ strain)

names(PSS_fact) <- row.names(PSS_factors)
ddsMat_PSS$sizeFactor <- PSS_fact
  
PSS_norm_counts <- counts(ddsMat_PSS, normalized=TRUE)
grp_File <- data.frame("WT"=apply(PSS_norm_counts[,c("WT1", "WT2", "WT3")],1,mean), "dWT"=apply(PSS_norm_counts[,c("dWT1", "dWT2", "dWT3")],1,mean), "TV"=apply(PSS_norm_counts[,c("TV1", "TV2")],1,mean))
 
grp_File$seqname <- rep(NA, length(grp_File$WT))
seqnames <- c("BA000022.2", "AP004310.1", "AP004311.1", "AP004312.1", "AP006585.1")
for(i in seqnames){
  grp_File[which(grepl(i, row.names(grp_File))),]$seqname <- i
}
grp_File$strand <- rep(NA, length(grp_File$WT))
for(i in c("plus", "minus")){
  grp_File[which(grepl(i, row.names(grp_File))),]$strand <- i
}

for(i in seqnames){
  for(j in c("plus", "minus")){
    tmp <- subset(grp_File, grp_File$seqname==i & grp_File$strand==j)[,c("WT", "dWT", "TV")]
    s=""
    if(j=="plus"){
      s="fw"
    }else{
      s="rev"
   }
    write.table(tmp, file=paste("Output/grp/PSS/", i, "_PSS-multireads_WT-dWT-TV_", s, ".grp", sep=""), sep="\t", row.names=FALSE, col.names = FALSE)
  }
}
  
rm(PSS_norm_counts)
}

1.2 TSS data

TSS_coverage <- read.table("Input/multireads_TSS_5ends_combined_5sensing.txt", sep="\t", header=TRUE, row.names=1)
ddsMat_TSS <- DESeqDataSetFromMatrix(countData = TSS_coverage,
                                 colData = coldata_general,
                                 design = ~ strain)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
TSS_fact <- TSS_factors$factor
names(TSS_fact) <- row.names(TSS_factors)
ddsMat_TSS$sizeFactor <- TSS_fact

TSS_norm_counts <- counts(ddsMat_TSS, normalized=TRUE)
rm(ddsMat_TSS)

if(FALSE){
grp_File <- data.frame("WT"=apply(TSS_norm_counts[,c("WT1", "WT2", "WT3")],1,mean), "dWT"=apply(TSS_norm_counts[,c("dWT1", "dWT2", "dWT3")],1,mean), "TV"=apply(TSS_norm_counts[,c("TV1", "TV2")],1,mean))
 
grp_File$seqname <- rep(NA, length(grp_File$WT))
seqnames <- c("BA000022.2", "AP004310.1", "AP004311.1", "AP004312.1", "AP006585.1")
for(i in seqnames){
  grp_File[which(grepl(i, row.names(grp_File))),]$seqname <- i
}
grp_File$strand <- rep(NA, length(grp_File$WT))
for(i in c("plus", "minus")){
  grp_File[which(grepl(i, row.names(grp_File))),]$strand <- i
}

for(i in seqnames){
  for(j in c("plus", "minus")){
    tmp <- subset(grp_File, grp_File$seqname==i & grp_File$strand==j)[,c("WT", "dWT", "TV")]
    s=""
    if(j=="plus"){
      s="fw"
    }else{
      s="rev"
   }
    write.table(tmp, file=paste("Output/grp/TSS/", i, "_TSS-multireads_WT-dWT-TV_", s, ".grp", sep=""), sep="\t", row.names=FALSE, col.names = FALSE)
  }
}
  
}

1.3 Transcript data

transcript_coverage <- read.table("Input/multireads_transcript_coverage_combined_5sensing.txt", sep="\t", header=TRUE, row.names=1)
ddsMat_transcript <- DESeqDataSetFromMatrix(countData = transcript_coverage,
                                 colData = coldata_general,
                                 design = ~ strain)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
transcript_fact <- transcript_factors$factor
names(transcript_fact) <- row.names(transcript_factors)
ddsMat_transcript$sizeFactor <- transcript_fact

transcript_norm_counts <- counts(ddsMat_transcript, normalized=TRUE)
rm(ddsMat_transcript)
rm(transcript_coverage)

if(FALSE){
  grp_File <- data.frame("WT"=apply(transcript_norm_counts[,c("WT1", "WT2", "WT3")],1,mean), "dWT"=apply(transcript_norm_counts[,c("dWT1", "dWT2", "dWT3")],1,mean), "TV"=apply(transcript_norm_counts[,c("TV1", "TV2")],1,mean))
  grp_File$seqname <- rep(NA, length(grp_File$WT))
  seqnames <- c("BA000022.2", "AP004310.1", "AP004311.1", "AP004312.1", "AP006585.1")
  for(i in seqnames){
    grp_File[which(grepl(i, row.names(grp_File))),]$seqname <- i
  }
  grp_File$strand <- rep(NA, length(grp_File$WT))
  for(i in c("plus", "minus")){
    grp_File[which(grepl(i, row.names(grp_File))),]$strand <- i
  }
  for(i in seqnames){
    for(j in c("plus", "minus")){
      tmp <- subset(grp_File, grp_File$seqname==i & grp_File$strand==j)[,c("WT", "dWT", "TV")]
      s=""
      if(j=="plus"){
        s="fw"
      }else{
        s="rev"
      }
      write.table(tmp, file=paste("Output/grp/transcript/", i, "_transcript-multireads_WT-dWT-TV_", s, ".grp", sep=""), sep="\t", row.names=FALSE, col.names = FALSE)
    }
  }
}

2 Perform DESeq2 analysis

names(PSS_coverage) <- paste(names(PSS_coverage), "_PSS", sep="")
names(PSS_fact) <- paste(names(PSS_fact), "_PSS", sep="")
names(TSS_coverage) <- paste(names(TSS_coverage), "_TSS", sep="")
names(TSS_fact) <- paste(names(TSS_fact), "_TSS", sep="")

# merge PSS and TSS to merged_raw
merged_raw <- merge(PSS_coverage, TSS_coverage, by="row.names")
row.names(merged_raw) <- merged_raw$Row.names
merged_raw$Row.names <- NULL
rm(TSS_coverage)

When taking into account lowly populated sites, DESeq2 dispersion estimates become over-fitted. Hence, edgeR::filterByExpression() is used as a pre-processing step before further DESeq2 analyses.

group=c(rep("dWT_PSS", 3), rep("WT_PSS", 3), rep("TV_PSS", 2), rep("dWT_TSS", 3), rep("WT_TSS", 3), rep("TV_TSS", 2))
y <- DGEList(counts=merged_raw, group=group)
nrow(y)
## [1] 7894038
keep <- filterByExpr(y)
y <- y[keep, ,keep.lib.size=FALSE]
nrow(y)
## [1] 19944
merged_filtered <- merged_raw[row.names(y$counts),]
nrow(merged_filtered)
## [1] 19944
rm(merged_raw)

coldata <- read.csv("Input/colData_PSS-TSS.csv", row.names=1)

# create DESeq2 data object
ddsMat_PSSTSS <- DESeqDataSetFromMatrix(countData = merged_filtered,
                                 colData = coldata,
                                 design = ~ strain + type)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
ddsMat_PSSTSS$sizeFactor <- c(PSS_fact, TSS_fact)
ddsMat_PSSTSS <- DESeq(ddsMat_PSSTSS)
## using pre-existing size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
res_PSSTSS <- results(ddsMat_PSSTSS, contrast=c("type", "PSS", "TSS"))
write.csv(res_PSSTSS[order(res_PSSTSS$padj),], file="Output/DESeq2_resultsTables/results_PSS-TSS-comparisons.csv")

2.1 Diagnostic Plots

plotDispEsts(ddsMat_PSSTSS)

pdf(file="Output/DESeq2_Plots/ddsMat_PSSTSS_DispEsts.pdf", width=4.5, height=4.5)
plotDispEsts(ddsMat_PSSTSS, xlab="Mean of Normalized Counts", ylab="Dispersion")
dev.off()
## png 
##   2
p <- PCA_plot_PSSTSS(ddsMat_PSSTSS, "PSS-TSS")
p

ggsave("Output/DESeq2_Plots/ddsMat_PSSTSS_PCA.pdf", plot=p, width=15, height=9, units="cm")

p <- heatmap_plot_PSSTSS(ddsMat_PSSTSS, "PSS-TSS")
p

ggsave("Output/DESeq2_Plots/ddsMat_PSSTSS_heatMap.pdf", plot=p, width=15, height=12, units="cm")
rm(ddsMat_PSSTSS) # to free space

pdf(file="Output/DESeq2_Plots/ddsMat_PSSTSS_hist_log2FC.pdf", width=4.5, height=4.5)
hist(res_PSSTSS$log2FoldChange, breaks=20, main="", xlab="Log2FC(PSS/TSS)", ylab="Frequency")
dev.off()
## png 
##   2
hist(res_PSSTSS$log2FoldChange, breaks=20, main="DESeq2 Normalization", xlab="Log2FC PSS/TSS", ylab="Frequency")

# non-normalized
log2_all <- log2(apply(merged_filtered[,1:8],1, mean)/apply(merged_filtered[,9:16],1,mean))
hist(log2_all, breaks=20, main="No normalization", xlab="Log2(PSS/TSS)")

rm(merged_filtered)

p <- MAplot_ggplot(res_PSSTSS, foldchange=0.8, y_axis_label = "Log2 fold-change(PSS/TSS)")
p <- p + scale_colour_manual(values=c("UP"="#4e9879ff", "DOWN"="#9b511fff", "NO"="#d3d3d3ff"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
p

ggsave("Output/DESeq2_Plots/ddsMat_PSSTSS_MAplot.pdf",plot=p, width=15, height=12, units="cm")

2.2 Extract DESeq2 results

res_PSSTSS <- subset(res_PSSTSS, !is.na(res_PSSTSS$padj))

# down: TSS, up: PSS
count_up_down(res_PSSTSS, foldchange=0.8, padjusted=0.05)
## [1] "number of features down:  11620"
## [1] "number of features up:  5424"
p <- volcanoPlot_ggplot(res_PSSTSS, foldchange=0.8, padjusted=0.05)
p

ggsave("Output/DESeq2_Plots/ddsMat_PSSTSS_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")

pdf(file="Output/DESeq2_Plots/ddsMat_PSSTSS_pValuePlot.pdf", width=4.5, height=4.5)
pvaluePlot(res_PSSTSS, "Comparison PSS to TSS")
dev.off()
## png 
##   2
pvaluePlot(res_PSSTSS, "Comparison PSS to TSS")

# extract which positions are TSS, which are PSS
PSS_positions_unfiltered <- row.names(subset(res_PSSTSS, res_PSSTSS$log2FoldChange>0.8 & res_PSSTSS$padj<0.05))
PSS_nonReduced_Ranges <- create_GRanges_object(base::strsplit(PSS_positions_unfiltered, "-"), res_PSSTSS[PSS_positions_unfiltered,]$baseMean)
TSS_positions_unfiltered <- row.names(subset(res_PSSTSS, res_PSSTSS$log2FoldChange<(-0.8) & res_PSSTSS$padj<0.05))

2.3 Filter for ratio 5’ ends / transcript coverage

trans_norm <- as.data.frame(transcript_norm_counts)
rm(transcript_norm_counts)
TSS_norm <- as.data.frame(TSS_norm_counts)
rm(TSS_norm_counts)
TSS_trans_ratio_cutoff <- 0.02
TSS_trans_ratio <- apply(TSS_norm[TSS_positions_unfiltered,],1,mean)/apply(trans_norm[TSS_positions_unfiltered,],1,mean)
TSS_above_cutoff <- subset(TSS_positions_unfiltered, TSS_trans_ratio[TSS_positions_unfiltered]>TSS_trans_ratio_cutoff)
TSS_positions_sorted<- base::strsplit(TSS_above_cutoff, "-")
TSS_nonReduced_Ranges <- create_GRanges_object(TSS_positions_sorted, res_PSSTSS[TSS_above_cutoff,]$baseMean)
rm(trans_norm)
rm(TSS_norm)

2.4 Create .gff files

PSS_positions_Ranges <- GenomicRanges::reduce(PSS_nonReduced_Ranges)
length(PSS_positions_Ranges)
## [1] 4821
TSS_positions_Ranges <- GenomicRanges::reduce(TSS_nonReduced_Ranges)
length(TSS_positions_Ranges)
## [1] 7406
TSS_positions_Ranges$type <- "TSS"
PSS_positions_Ranges$type <- "PSS"
save_gff(c(TSS_positions_Ranges, PSS_positions_Ranges), "Output/gffs/", "TSS_PSS_multireads")
rm(TSS_positions_Ranges)
rm(PSS_positions_Ranges)

3 Perform DESeq2 analysis for PSS

PSS_use <- PSS_coverage[PSS_positions_unfiltered,]
rm(PSS_coverage)

coldata <- read.csv("Input/colData_PSS.csv", row.names=1)

ddsMat_PSS <- DESeqDataSetFromMatrix(countData = PSS_use,
                                 colData = coldata,
                                 design = ~ strain)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
ddsMat_PSS$sizeFactor <- PSS_fact
ddsMat_PSS <- DESeq(ddsMat_PSS) 
## using pre-existing size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
rm(PSS_use)

3.1 Diagnostic Plots

plotDispEsts(ddsMat_PSS, main="PSS comparison", xlab="Mean of Normalized Counts", ylab="Dispersion")

p <- PCA_plot(ddsMat_PSS, "PSS")
p 

ggsave("Output/DESeq2_Plots/ddsMat_PSS_PCA.pdf", plot=p, width=9, height=9, units="cm")
p <- heatmap_plot(ddsMat_PSS, "PSS")
p

ggsave("Output/DESeq2_Plots/ddsMat_PSS_heatMap.pdf", plot=p, width=15, height=12, units="cm")

3.2 Extract Results

# extract results
PSS_result_dWT_WT <- results(ddsMat_PSS, contrast=c("strain", "dWT", "WT")) # dWT/WT -> higher in dWT: higher log2FC
write.csv(PSS_result_dWT_WT[order(PSS_result_dWT_WT$padj),], file="Output/DESeq2_resultsTables/results_PSS-dWT-WT.csv")

PSS_result_dWT_TV <- results(ddsMat_PSS, contrast=c("strain", "dWT", "TV")) # dWT/TV -> higher in dWT: higher log2FC
write.csv(PSS_result_dWT_TV[order(PSS_result_dWT_TV$padj),], file="Output/DESeq2_resultsTables/results_PSS-dWT-TV.csv")

PSS_result_WT_TV <- results(ddsMat_PSS, contrast=c("strain", "WT", "TV")) # WT/TV -> higher in WT: higher log2FC
write.csv(PSS_result_WT_TV[order(PSS_result_WT_TV$padj),], file="Output/DESeq2_resultsTables/results_PSS-WT-TV.csv")

rm(ddsMat_PSS)
indices_plus <- which(strand(PSS_nonReduced_Ranges)=="+")
indices_minus <- which(strand(PSS_nonReduced_Ranges)=="-") 
PSS_names <- c()
PSS_names[indices_plus] <- paste(seqnames(PSS_nonReduced_Ranges[indices_plus]),start(PSS_nonReduced_Ranges[indices_plus]),"plus",sep="-")
PSS_names[indices_minus] <- paste(seqnames(PSS_nonReduced_Ranges[indices_minus]),start(PSS_nonReduced_Ranges[indices_minus]),"minus",sep="-")
rm(indices_plus)
rm(indices_minus)
PSS_nonReduced_Ranges$names <- PSS_names
PSS_nonReduced_Ranges$featuresOverlap <- rep("", length(PSS_nonReduced_Ranges))
PSS_nonReduced_Ranges$TU_overlap <- rep("", length(PSS_nonReduced_Ranges))
for(i in 1:length(PSS_nonReduced_Ranges)){
  PSS_nonReduced_Ranges$featuresOverlap[i] <- paste(features[which(countOverlaps(features,PSS_nonReduced_Ranges[i])>0)]$locus_tag,collapse=",")
  PSS_nonReduced_Ranges$TU_overlap[i] <- paste(TUs[which(countOverlaps(TUs,PSS_nonReduced_Ranges[i])>0)]$index,collapse=",")
}
# create tables with annotation: which features and TUs are overlapping with PSS?
df_PSS_nonRed_Ranges <- as.data.frame(PSS_nonReduced_Ranges)
row.names(df_PSS_nonRed_Ranges) <- df_PSS_nonRed_Ranges$names

# PSS result dWT WT
PSS_result_dWT_WT_annot <- rownames_to_column(as.data.frame(PSS_result_dWT_WT[order(PSS_result_dWT_WT$padj),]))
PSS_result_dWT_WT_annot$CDS <- df_PSS_nonRed_Ranges[PSS_result_dWT_WT_annot$rowname,]$featuresOverlap
PSS_result_dWT_WT_annot$TUs <- df_PSS_nonRed_Ranges[PSS_result_dWT_WT_annot$rowname,]$TU_overlap
PSS_result_dWT_WT_annot <- PSS_result_dWT_WT_annot[,c(8,9,1,2:7)]
write_tsv(PSS_result_dWT_WT_annot, "Output/DESeq2_resultsTables/results_PSS_dWT_WT_annotated.tsv")

# PSS result dWT TV
PSS_result_dWT_TV_annot <- rownames_to_column(as.data.frame(PSS_result_dWT_TV[order(PSS_result_dWT_TV$padj),]))
PSS_result_dWT_TV_annot$CDS <- df_PSS_nonRed_Ranges[PSS_result_dWT_TV_annot$rowname,]$featuresOverlap
PSS_result_dWT_TV_annot$TUs <- df_PSS_nonRed_Ranges[PSS_result_dWT_TV_annot$rowname,]$TU_overlap
PSS_result_dWT_TV_annot <- PSS_result_dWT_TV_annot[,c(8,9,1,2:7)]
write_tsv(PSS_result_dWT_TV_annot, "Output/DESeq2_resultsTables/results_PSS_dWT-TV_annotated.tsv")

# PSS result WT, TV
PSS_result_WT_TV_annot <- rownames_to_column(as.data.frame(PSS_result_WT_TV[order(PSS_result_WT_TV$padj),]))
PSS_result_WT_TV_annot$CDS <- df_PSS_nonRed_Ranges[PSS_result_WT_TV_annot$rowname,]$featuresOverlap
PSS_result_WT_TV_annot$TUs <- df_PSS_nonRed_Ranges[PSS_result_WT_TV_annot$rowname,]$TU_overlap
PSS_result_WT_TV_annot <- PSS_result_WT_TV_annot[,c(8,9,1,2:7)]
write_tsv(PSS_result_WT_TV_annot, "Output/DESeq2_resultsTables/results_PSS_WT-TV_annotated.tsv")

Run in respective Directory (output/DESeq2_resultsTables/) python changeAnnotation_DESeq2.py results_PSS_dWT_WT_annotated.tsv results_dWT_WT_annotated-2.tsv python changeAnnotation_DESeq2.py results_PSS_dWT-TV_annotated.tsv results_dWT-TV_annotated-2.tsv python changeAnnotation_DESeq2.py results_PSS_-_WT-TV_annotated.tsv results_PSS-WT-TV_annotated-2.tsv

pdf(file="Output/DESeq2_Plots/ddsMat_PSS-dWT-WT_pValuePlot.pdf", width=4.5, height=4.5)
pvaluePlot(PSS_result_dWT_WT, "PSS dWT WT")
dev.off()
## png 
##   2
pdf(file="Output/DESeq2_Plots/ddsMat_PSS-dWT-TV_pValuePlot.pdf", width=4.5, height=4.5)
pvaluePlot(PSS_result_dWT_TV, "PSS dWT dTV")
dev.off()
## png 
##   2
pvaluePlot(PSS_result_dWT_WT, "PSS dWT WT")

pvaluePlot(PSS_result_dWT_TV, "PSS dWT dTV")

Filter out PSS below filter level of results(). For these, p.adj is set to NA - hence: remove positions at which p.adj==NA

PSS_result_dWT_WT <- subset(PSS_result_dWT_WT, !is.na(PSS_result_dWT_WT$padj))
PSS_result_dWT_TV <- subset(PSS_result_dWT_TV, !is.na(PSS_result_dWT_TV$padj))
PSS_result_WT_TV <- subset(PSS_result_WT_TV, !is.na(PSS_result_WT_TV$padj))

3.2.1 PSS dWT WT

count_up_down(PSS_result_dWT_WT, foldchange=1, padjusted=0.05)
## [1] "number of features down:  4"
## [1] "number of features up:  67"
p <- volcanoPlot_ggplot(as.data.frame(PSS_result_dWT_WT), foldchange=1, padjusted=0.05) + scale_colour_manual(values=c("DOWN"="#000000ff", "UP"="#005a96ff", "NO"="#d3d3d3ff"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
p

ggsave("Output/DESeq2_Plots/ddsMat_PSS-0h_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")

p <- MAplot_ggplot(PSS_result_dWT_WT, foldchange=1.0, y_axis_label = "Log2 fold-change(rne(WT)/WT)") + scale_colour_manual(values=c("DOWN"="#000000ff", "UP"="#005a96ff", "NO"="#d3d3d3ff"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
p

ggsave("Output/DESeq2_Plots/ddsMat_PSS-dWT-WT_MAplot.pdf",plot=p, width=15, height=12, units="cm")

3.2.2 PSS dWT-TV

count_up_down(PSS_result_dWT_TV, foldchange=1, padjusted=0.05)
## [1] "number of features down:  175"
## [1] "number of features up:  94"
p <- volcanoPlot_ggplot(as.data.frame(PSS_result_dWT_TV), foldchange=1, padjusted=0.05) + scale_colour_manual(values=c("DOWN"="#e69f00ff", "UP"="#005a96ff", "NO"="#d3d3d3ff"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
p

ggsave("Output/DESeq2_Plots/ddsMat_PSS-dWT-TV_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")

p <- MAplot_ggplot(PSS_result_dWT_TV, foldchange=1.0, y_axis_label = "Log2 fold-change(rne(WT)/rne(5p))") + scale_colour_manual(values=c("DOWN"="#e69f00ff", "UP"="#005a96ff", "NO"="#d3d3d3ff"))# + ylim(-5,+5)
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
p

ggsave("Output/DESeq2_Plots/ddsMat_PSS-1h_MAplot.pdf",plot=p, width=15, height=12, units="cm")

3.2.3 PSS WT-TV

count_up_down(PSS_result_WT_TV, foldchange = 1, padjusted = 0.05)
## [1] "number of features down:  237"
## [1] "number of features up:  111"
p <- MAplot_ggplot(PSS_result_WT_TV, foldchange=1.0, y_axis_label = "Log2 fold-change(WT/rne(5p))") + scale_colour_manual(values=c("DOWN"="#e69f00ff", "UP"="#000000ff", "NO"="#d3d3d3ff"))# + ylim(-5,+5)
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
p

ggsave("Output/DESeq2_Plots/ddsMat_PSS-WT-TV_MAplot.pdf",plot=p, width=15, height=12, units="cm")

Session info

## R version 4.0.4 (2021-02-15)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 11 (bullseye)
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
## 
## locale:
##  [1] LC_CTYPE=de_DE.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=de_DE.UTF-8        LC_COLLATE=de_DE.UTF-8    
##  [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=de_DE.UTF-8   
##  [7] LC_PAPER=de_DE.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] viridis_0.5.1               viridisLite_0.3.0          
##  [3] rtracklayer_1.50.0          Biostrings_2.58.0          
##  [5] XVector_0.30.0              vsn_3.58.0                 
##  [7] RColorBrewer_1.1-2          pheatmap_1.0.12            
##  [9] ggpubr_0.4.0                ggrepel_0.9.1              
## [11] edgeR_3.32.1                limma_3.46.0               
## [13] DESeq2_1.30.1               SummarizedExperiment_1.20.0
## [15] Biobase_2.50.0              MatrixGenerics_1.2.1       
## [17] matrixStats_0.58.0          GenomicRanges_1.42.0       
## [19] GenomeInfoDb_1.26.2         IRanges_2.24.1             
## [21] S4Vectors_0.28.1            BiocGenerics_0.36.0        
## [23] forcats_0.5.1               stringr_1.4.0              
## [25] dplyr_1.0.5                 purrr_0.3.4                
## [27] readr_1.4.0                 tidyr_1.1.3                
## [29] tibble_3.1.0                ggplot2_3.3.3              
## [31] tidyverse_1.3.0             knitr_1.31                 
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_2.0-0         ggsignif_0.6.1           ellipsis_0.3.1          
##  [4] rio_0.5.26               fs_1.5.0                 rstudioapi_0.13         
##  [7] farver_2.1.0             affyio_1.60.0            bit64_4.0.5             
## [10] AnnotationDbi_1.52.0     fansi_0.4.2              lubridate_1.7.10        
## [13] xml2_1.3.2               splines_4.0.4            cachem_1.0.4            
## [16] geneplotter_1.68.0       jsonlite_1.7.2           Rsamtools_2.6.0         
## [19] broom_0.7.5              annotate_1.68.0          dbplyr_2.1.0            
## [22] BiocManager_1.30.10      compiler_4.0.4           httr_1.4.2              
## [25] backports_1.2.1          assertthat_0.2.1         Matrix_1.3-2            
## [28] fastmap_1.1.0            cli_2.3.1                htmltools_0.5.1.1       
## [31] tools_4.0.4              affy_1.68.0              gtable_0.3.0            
## [34] glue_1.4.2               GenomeInfoDbData_1.2.4   Rcpp_1.0.6              
## [37] carData_3.0-4            cellranger_1.1.0         jquerylib_0.1.3         
## [40] vctrs_0.3.6              preprocessCore_1.52.1    xfun_0.22               
## [43] openxlsx_4.2.3           rvest_1.0.0              lifecycle_1.0.0         
## [46] rstatix_0.7.0            XML_3.99-0.5             zlibbioc_1.36.0         
## [49] scales_1.1.1             hms_1.0.0                yaml_2.2.1              
## [52] curl_4.3                 gridExtra_2.3            memoise_2.0.0           
## [55] sass_0.3.1               stringi_1.5.3            RSQLite_2.2.3           
## [58] highr_0.8                genefilter_1.72.1        zip_2.1.1               
## [61] BiocParallel_1.24.1      rlang_0.4.10             pkgconfig_2.0.3         
## [64] bitops_1.0-6             evaluate_0.14            lattice_0.20-41         
## [67] labeling_0.4.2           GenomicAlignments_1.26.0 bit_4.0.4               
## [70] tidyselect_1.1.0         magrittr_2.0.1           R6_2.5.0                
## [73] generics_0.1.0           DelayedArray_0.16.2      DBI_1.1.1               
## [76] pillar_1.5.1             haven_2.3.1              foreign_0.8-81          
## [79] withr_2.4.1              abind_1.4-5              survival_3.2-7          
## [82] RCurl_1.98-1.2           modelr_0.1.8             crayon_1.4.1            
## [85] car_3.0-10               utf8_1.1.4               rmarkdown_2.7           
## [88] locfit_1.5-9.4           grid_4.0.4               readxl_1.3.1            
## [91] data.table_1.14.0        blob_1.2.1               reprex_1.0.0            
## [94] digest_0.6.27            xtable_1.8-4             munsell_0.5.0           
## [97] bslib_0.2.4